* CODE TO RUN SIMULATIONS/BUILD SIM DATA SET


********************************************************************************
** Set Simulation Parameters
********************************************************************************

* Specify form of tax vector
local uniformTax 1

* Specify value of tax (if uniform) in $ per ton
local carbonTax 25

* Specify foreign emissions intensity (1= foreign; 2=domestic proxy)
local foreign_ei 2

forvalues tax=25(1)45 {
foreach scenario in "tax" "subsidy" "subsidy_h_p85" "subsidy_e_p90" {

********************************************************************************
** Initialize the results file
********************************************************************************
local filename = "`scenario'"

*					 
********************************************************************************
** I. Read in the elasticities for each specifications
********************************************************************************	
{
			  
	* 1. import regression results
		
	use "$result_regs/regression_results_wide.dta", clear	
	
	* renaming
	rename betaS beta_linei_prod
	rename betaM beta_linei_imp
	rename betaX beta_linei_exp
	
	* renaming ei2
	rename beta2S beta_linei2_prod
	rename beta2M beta_linei2_imp
	rename beta2X beta_linei2_exp
		
	
	* labels
	label var beta_linei_prod "linear coeff for production"
	label var beta_linei_imp "linear coeff for imports"
	label var beta_linei_exp "linear coeff for exports"

	label var beta_linei2_prod "quadratic coeff for production"
	label var beta_linei2_imp "quadratic coeff for imports"
	label var beta_linei2_exp "quadratic coeff for exports"
	
	keep year spec eivar model beta_linei*
}

********************************************************************************
** II. Import simulation data and merge/create elasticities & transfer rates
********************************************************************************
{
	joinby using "$result_simul/simulation_dataset.dta"

	
* *****************    generate domestic energy price elasticities

	gen beta_prod = beta_linei_prod * ei
	gen beta_imp = beta_linei_imp * ei
	gen beta_exp = beta_linei_exp * ei
	
	replace beta_prod = beta_linei_prod * ei_tot if eivar == "ei_tot"
	replace beta_imp = beta_linei_imp * ei_tot if eivar == "ei_tot"
	replace beta_exp = beta_linei_exp * ei_tot if eivar == "ei_tot"
	
	replace beta_prod = beta_prod +  beta_linei2_prod * (ei^2) if eivar=="ei" & model=="lin_ei2"
	replace beta_imp = beta_imp +  beta_linei2_imp * (ei^2) if eivar=="ei" & model=="lin_ei2"
	replace beta_exp = beta_exp +  beta_linei2_exp * (ei^2) if eivar=="ei" & model=="lin_ei2"

	replace beta_prod = beta_prod +  beta_linei2_prod * (ei_tot^2) if eivar=="ei_tot" & model=="lin_ei2"
	replace beta_imp = beta_imp +  beta_linei2_imp * (ei_tot^2) if eivar=="ei_tot" & model=="lin_ei2"
	replace beta_exp = beta_exp +  beta_linei2_exp * (ei_tot^2) if eivar=="ei_tot" & model=="lin_ei2"
	
	label var beta_prod "elasticity for production"
	label var beta_imp "elasticity for imports"
	label var beta_exp "elasticity for exports"



********************   merge in transfer rate   

	merge 1:1 naics spec using  "$result_regs/regression_results_tr_by_ind.dta", nogen keep(1 3)

********************  set carbon tax ***************************
	
	* carbon tax
	gen carbonTax = `tax'
	label var carbonTax "$ per ton"

	duplicates drop
}


********************************************************************************
** ****************III. Build price changes implied by taxes and subsidies
**************************************************************************************
{
	
******************** 1. Calculate implied industry-level price increase for given carbon tax on metric tons/mmbtu 
	
	*(1) caculate price change with a tax
	if "`scenario'" == "tax" {
		gen change_p = (ci_tvshare / 1000) * carbonTax // CarbonI_TV units=kg/mmbtu, carbon tax in $ per ton
		gen subsidy_rate = 0.0
		gen subsidy = 0.0
	}
	
	*(2) caculate subsidy under optimal rates and price change	
	if "`scenario'" != "tax" {
	
		* (a) calculate subsidy rate
		if "`scenario'" == "subsidy" {
			egen subsidy_rate = median(transfer_rate), by(naics eivar)
		}
		if "`scenario'" == "subsidy_h" {
			egen subsidy_rate = median(tr_qreg), by(naics eivar)
		}
		if (strmatch("`scenario'","subsidy*_p??")) {
			
			* determine elegibility threshold based on leakage rate
			* define threshold
			if (strmatch("`scenario'","*_e_*")) {
				summ e_d, det
				local limit_l = r(p25)
				local limit_u = r(p90)
				summ te_2007, det
				local limit_te = r(p90)
				gen elegible = 0
				replace elegible = 1 if (te_2007 > `limit_te' & e_d > `limit_l') |  (e_d > `limit_u')
			}
			else {
				if (strmatch("`scenario'","*_h*") | strmatch("`scenario'","*_ih*")) {
					egen median_lr = median(lr_qreg), by(naics eivar)
				}
				else if (strmatch("`scenario'","*_i*")) {
					egen median_lr = median(leakage_rate), by(naics eivar)
				}
				else {
					egen median_lr = median(leakage_rate), by(naics eivar)
				}				
				local limit = substr("`scenario'",-2,2)
				_pctile median_lr if eivar=="ei", p(`limit')
				gen elegible = 1 if median_lr >= r(r1) & eivar=="ei"
				replace elegible = 0 if median_lr < r(r1) & eivar=="ei"
				_pctile median_lr if eivar=="ei_tot", p(`limit')
				replace elegible = 1 if median_lr >= r(r1) & eivar=="ei_tot"
				replace elegible = 0 if median_lr < r(r1) & eivar=="ei_tot"
				drop median_lr
			}			
			
			* define subsidy rate
			if (strmatch("`scenario'","*_i_*")) {
				egen subsidy_rate = median(transfer_rate), by(naics eivar)
			}
			else if (strmatch("`scenario'","*_ih_*")) {
				egen subsidy_rate = median(tr_qreg), by(naics eivar)
			}
			else if (strmatch("`scenario'","*_h_*")) {
				* egen subsidy_rate = median(tr_qreg) if elegible==1, by(eivar)
				gen subsidy_rate = .80 if elegible==1 
			}
			else if (strmatch("`scenario'","*_e_*")) {
				* egen subsidy_rate = median(tr_qreg) if elegible==1, by(eivar)
				gen subsidy_rate = .80 if elegible==1 
			}
			else {
				egen subsidy_rate = median(transfer_rate) if elegible==1, by(eivar)
			}
			replace subsidy_rate = 0 if elegible==0
		}
		
		* (b) calculate subsidy
		gen subsidy = subsidy_rate * 25 * e_d if eivar=="ei"
		replace subsidy= subsidy_rate * 25 * e_f if eivar=="ei_tot"

		* (c) calculate change in price
		gen change_p = ci_tvshare / 1000 * carbonTax 
		replace change_p = change_p - (subsidy / ei_tv)
		
	}
	
******************* now construct price changes for the comprehensive energy prices.	
	{
	
	* if ei_tot specification, convert to indirect price change
	tempfile indirect_p_change
	
	* (1) get naics2007-2012 crosswalks & naics2007-specid panel
	preserve
		do "$analyzepath/code/subfunctions/crosswalk_naics2007-2012.do"
		tempfile crosswalk
		save `crosswalk.dta', replace
	restore
	
	* (2) convert
		* get price change
	preserve
		keep naics change_p shipments spec
		rename naics naics2012
		rename change_p change_p
		duplicates drop
		tempfile change_p
		save `change_p.dta', replace
	restore
	
	preserve
	
		* merge in naics2007 code
		sum spec
		local N = `r(max)'
		use `crosswalk.dta', clear
		expand `N'
		bysort naics2007: gen spec = _n
		sort naics2012 naics2007 spec
		merge m:1 naics2012 spec using `change_p.dta', keep(2 3) nogen
		
		* convert to change_p_ind 
		rename naics2007 naics
		do "$analyzepath/code/subfunctions/convert_indirectprice.do"
		
		* bring back naics2012
		rename naics naics2007
		merge m:1 naics2007 using `crosswalk.dta', keep(1 3) nogen
		
		replace change_p_ind = change_p_ind*shipments
		collapse (sum) change_p_ind shipments (first) change_p, by(naics2012 spec)
		replace change_p_ind = change_p_ind/shipments
		drop shipments
		
		rename naics2012 naics
		save `indirect_p_change', replace
	restore
	
	* (3) merge
	
	drop change_p
	
	drop if spec==.
	merge m:1 naics spec using `indirect_p_change', nogen
	
	gen pct_change_p = change_p / price_TIV
	replace pct_change_p = change_p_ind / price_TIV if eivar == "ei_tot"
	
	}
}

********************************************************************************
** ****************IV. Summarize changes
**************************************************************************************

** Each observation corresponds to a specification and naics.
** For each specification, we can compute prodn changes, trade changes, abatement and leakage.
** Following code computes outcome at the level of specification
qui {

	gen nettrade=imports + shipments
		
	* *******************2. Production and consumption changes	(subsidy specific)
	
	gen pct_change_prod = beta_prod * pct_change_p
	gen change_prod = pct_change_prod * shipments 
	*this is DP_i... units = millions of USD?

	********************* 3) Generate changes in foreign production

	* import/export changes
	
	gen pct_change_imp = beta_imp * pct_change_p 
	gen pct_change_exp = beta_exp * pct_change_p 
	
	
	gen change_imp = pct_change_imp * imports
	gen change_exp = pct_change_exp * exports
		
		
	********************* 4) calculate transfer rate (check) and summarize
		
	gen trate = (change_imp - change_exp) / (-change_prod)
	label var trate"median of transfer rate check (first take median of elasticities, and then calculate transfer rate)"
  
		
	***********************5) Abatement

	gen abatement = change_prod * e_d 
	label var abatement "industry specific abatement in tons"
	
	gen abaterate =abatement/emissions_bau
	label var abaterate "industry specific abatement percentage"

	***************************6. transfer rate summaries
	*** lots of ways to summarize the range of transfer rate estimates across specifications nd industries	
	foreach tr in "tr_reg" "tr_qreg" "tr_p50" "tr_avg" {	
		bysort spec: egen `tr'_spec_mn=mean(`tr')     
	}
	
	*******************7. leakage 
	
	gen leakage = change_imp * e_d - change_exp * e_d
	replace leakage = change_imp * e_f - change_exp * e_f if eivar=="ei_tot"
	label var leakage "industry specific emission leakage in tons"
	
	
	gen lr= leakage/abatement
		
	************************8. net abatement
	gen netabate= (abatement+ leakage)/1000000.0
	
	gen total_emi_bau= emissions_bau + imports* e_d - exports * e_d
	replace total_emi_bau= emissions_bau + imports* e_f - exports * e_f if eivar=="ei_tot"

	gen change_emi = change_prod * e_d - (change_imp * e_d) + (change_exp * e_d)
	replace  change_emi = change_prod * e_d - (change_imp * e_f) + (change_exp * e_f) if eivar=="ei_tot"
	
	*************7. tax  revenues used as subsidies 	
	gen tax_rev=carbonTax  * e_d * (shipments + change_prod) 	
	* million USD

	label var tax_rev "industry specific tax revenues, $m"
	gen subsidy_paid= subsidy * (shipments + change_prod) 

	
	
	gen consumption=shipments-exports+imports

	 *****************
	 
	* Now Compute aggregate values - either summed w/i specification or summarized w/i specification
	bysort spec: egen leakage_agg 			= total(leakage/1000000.0)
	label var leakage_agg 		"aggregate leakage in M tons"
	
	
	bysort spec: egen emission_bau_agg		= total(emissions_bau/1000000.0)
	label var emission_bau_agg 	"aggregate BAU emission in M tons "
	
	
	bysort spec: egen abatement_agg 		= total(abatement/1000000.0)
	label var abatement_agg 	"aggregate emission abatement in M tons "
	
	bysort spec: egen totemi_bau_agg		= total(total_emi_bau/1000000.0)
	label var totemi_bau_agg 	"aggregate BAU emission (global) in M tons "
	
	
	bysort spec: gen leakagerate_agg 		=  - leakage_agg/abatement_agg 
	label var leakagerate_agg 	"aggregate leakage rate"
		

	bysort spec: gen abatementrate_agg		= abatement_agg/emission_bau_agg
	label var abatementrate_agg "aggregate abatement rate"
		
	
	bysort spec: egen change_prod_agg 		= total(change_prod/1000.0)
	label var change_prod_agg 	"change in shipments in $1m, gdpdef"
	

	bysort spec: egen change_imp_agg 		= total(change_imp/1000.0)
	label var change_imp_agg 	"change in imports in $1m, gdpdef"
	
	bysort spec: egen change_exp_agg 		= total(change_exp/1000.0)
	label var change_exp_agg 	"change in exports in $1m, gdpdef"
	
	bysort spec: egen prod_agg 				= total(shipments/1000.0)
	label var prod_agg 			"value of shipments in $1b, gdpdef"
	
	bysort spec: egen imp_agg 				= total(imports/1000.0)
	label var imp_agg 			"value of imports in $1b, gdpdef"
	
	bysort spec: egen exp_agg 				= total(exports/1000.0)
	label var exp_agg 			"value of exports in $1b, gdpdef"
	
	bysort spec: egen net_agg 				= total(nettrade/1000.0)
	label var net_agg 			"value of net exports in $1b, gdpdef"
	
	bysort spec: egen tr_agg = median(subsidy_rate) if subsidy!=0
	sort spec tr_agg
	by spec: replace tr_agg = tr_agg[1] if tr_agg==.
	replace tr_agg = 0 if tr_agg==.
	label var tr_agg 	"median transfer rate conditional on subsidy"
	
	bysort spec: egen tax_rev_agg 			= total(tax_rev)
	label var tax_rev_agg 	"aggregate tax revenues, $m"
	
	
	bysort spec: gen pc_change_prod_agg 		=  change_prod_agg / prod_agg
	
	bysort spec: gen pc_change_imp_agg 		=  change_imp_agg / imp_agg
	
	bysort spec: gen pc_change_exp_agg 		=  change_exp_agg / exp_agg
	
	bysort spec: egen subsidy_agg           = total(subsidy_paid) 
	label var subsidy_agg "total subsidy paid"
	
	bysort spec: egen netabate_agg           = total(netabate)
	
	
	bysort spec: gen netabatetrate_agg		= netabate_agg/totemi_bau_agg	
	label var abatementrate_agg "aggregate abatement rate"
	
	
	bysort spec: gen netabate_agg_check      = abatement_agg + leakage_agg 
	
	bysort spec: egen tr_reg_average = mean(tr_avg)
	
	bysort spec: gen net_tax_rev_      = tax_rev_agg - subsidy_agg  


}
*
	 
* save	
save "$result_simul/simulation_allspec_`filename'_`tax'.dta", replace
clear

}
}

